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Abstract. The distribution of partition function zeros is studied for the ± J model 
of spin glasses on the Bethe lattice. We find a relation between the distribution of 
complex cavity fields and the density of zeros, which enables us to obtain the density 
of zeros for the infinite system size by using the cavity method. The phase boundaries 
thus derived from the location of the zeros are consistent with the results of direct 
analytical calculations. This is the first example in which the spin-glass transition is 
related to the distribution of zeros directly in the thermodynamical limit. We clarify 
how the spin-glass transition is characterized by the zeros of the partition function. 
It is also shown that in the spin-glass phase a continuous distribution of singularities 
touches the axes of real field and temperature. 
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1. Introduction 

The physical quantities of a spin system arc completely determined by the location of 
the zeros of partition function on the complex parameter planes, such as the temperature 
and external field. This means that all the singularities of the system should be identified 
with zeros, and a phase transition occurs when a system parameter is driven across an 
accumulation point of zeros. 

For ferromagnetic systems, Lee and Yang proved that the zeros on the complex 
H plane (Lee- Yang zeros) lie on the imaginary axis. This result is the famous circle 
theorem [1]. A ferromagnetic transition occurs when the zeros touch the real-field axis 
aX H = 0, and zeros split the complex-field plane in two regions, > and < 0, 
where 3? denotes the real part. The density of zeros at the origin is directly proportional 
to the spontaneous magnetization in the ferromagnetic phase. As for the temperature 
T plane (Fisher zeros). Fisher found that zeros of a pure ferromagnetic system on the 
square lattice in the absence of an external field lie on the unit circle in the complex 
plane of sinh(2^J), where J > is the couphng constant and P — 1/T the inverse 
temperature [2]. 

The Lcc-Yang circle theorem is valid for random ferromagnets (where all the 
interactions Jij arc greater than or equal to 0) and all the Lcc-Yang zeros lie on 
the imaginary-field axis. Therefore the problem is reduced to finding the density 
of zeros on the imaginary- field axis g{0), which can be calculated from the analytic 
continuation of the magnetization as a function of the real positive field m{H), as 
27rg{9) — 3? m{2/3H :— iO) (— tt < ^ < tt) [1]. For diluted ferromagnets, several works 
have explored the connection between the Griffiths singularity [3] and the density of 
zeros both in solvable models [4, 5] and experiments [6]. These studies suggest that the 
density of zeros has an essential singularity on the imaginary field axis at the origin as 
g{0) ~ exp (— A/I^^l), where A is a positive constant. The essential singularity at the 
origin, in Griffiths' interpretation [3] , is caused by the presence of large clusters due to 
random fiuctuations of the interactions. 

For spin glasses, the situation is more complicated. The locations of zeros are 
not restricted to the imaginary-field axis. Since it is difficult to treat analytically 
the partition functions of these systems, zeros have been studied mainly by numerical 
evaluation of partition functions of finite-size systems on the complex external field [7, 8] 
and temperature [9, 10, 11] planes. The zeros were found to be distributed generally 
on two-dimensional areas, not only along a line as in ferromagnets. Recently a detailed 
study of the density of zeros on the imaginary-field axis suggested the existence of the 
Griffiths singularity also in spin-glass models in the paramagnetic phase [8]. 

It is therefore highly desirable to study the distribution of zeros of spin glasses in the 
limit of infinite system size since the complex behaviour of spin glasses should manifest 
itself in non-trivial distributions of zeros. We therefore investigate in this paper the 
relation between the distribution of zeros and the spin-glass transition in the infinite 
system size for the spin-glass model on the Bethe lattice. We here define the Bethe 
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lattice as the interior part of the infinitely large Cayley tree. This definition enables us 
to investigate the zeros of infinite size systems without approximations (even within the 
spin-glass phase) using the cavity method [12], which is perfectly suited for our purpose 
of studying the distribution of zeros of spin glasses. From the distribution of zeros we 
successfully detect the phase boundary of the spin-glass phase, defined as the line where 
the spin-glass susceptibility diverges, in accordance with previous studies. We also show 
that the system is singular everywhere as a function of the real- valued temperature or 
field in the spin-glass phase. 

The outhne of this paper is as follows. We first introduce the models and the 
relationship between the density of zeros and the distribution of cavity field in the next 
section. Then we explain the actual numerical procedures in section 3. The results are 
described in section 4, where distributions of zeros on the complex field and temperature 
planes are investigated and the phase diagrams are determined. The existence of the 
Griffiths singularity is also discussed. Finally we conclude this paper in section 5. 

2. Formulation 

In this section, the main ideas of this paper are presented. It is shown that the zeros 
of the partition function for the ± J model on the Bethe lattice are efficiently evaluated 
using the cavity method. 

2.1. Cayley tree and the cavity method 

Let us consider the ± J model of spin glasses on a Cayley tree in a uniform magnetic field 
H. The Cayley tree is a cycle- free graph where each site is connected to c neighbours 
(figure 1). The Hamiltonian is 



where {i,j) is a nearest neighbour pair, and the value of interaction Jij distributes as 




(1) 



P (Jij) = p5 {Jij -J) + {l-p)S {Jij + J) 



(2) 



with < p < 1 and J — 1. 




Figure 1. The local structure of the tree system with the coordination number c = 3. 
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On the Cayley tree we can efficiently calculate the partition function by iteratively 
summing over spin variables layer by layer from the outer boundary. The partial 
partition function Zi, which represents the partition function in the absence of deeper 
layers than the site i, is updated as 

c-1 

Zi = 2 cosh (/3/ii)n 

where j labels c — 1 spins on the layer previous to the current site i (figure 1). The 
independent effects of c — 1 spins on the previous layer have been passed to the present 
layer in terms of the cavity field hi and cavity biases {uj}, the definitions of which are 
given by 




c-1 

^H + Y.y'jiJij^hj). (5) 

i=i 

We hereafter assume that the function tanh^^ takes the principal branch, which restricts 
the value of the imaginary part of the bias Uj to a range [— 7r//3, n/fS]. From now on we 
assume the imaginary part of all the fields defined modulo 7r//3 and the principal branch 
of the tanh~^ is considered. The reader should keep this in mind, as we will not burden 
the notation with this condition explicitly. 

At the final step of this procedure, we consider the contribution from the central 
site of the Cayley tree, and the partition function of the whole system Z^'^^ is obtained 
as 

c 

= 2cosh (/3/i^'^)n 

and 

c 

h^^^ = H + J2uAJ.„h,), (7) 

where the superscript ^'^^ represents the central site. 

To consider the typical behaviour of the Cayley tree (where we assume uncorrelated, 
typical boundary conditions, to be specified explicitly below), we perform the average 
over the quenched randomness, which introduces the distribution of the cavity field at 
the Ith layer -P^«. The update rule of Pjjg is given by 




where [■ ■ ■] j denotes the average over the random interactions J^. In the thermodynamic 
limit, and if the distributions g converge as / ^ oo to PH,j3 = lim;_^oo Pr bi the 



cosh (/3 J, 



13) 



cosh ifiuj ( Jij ,hj)) 



Z, 



(3) 



cosh {13 Jij) 
cosh (f3uj(Jij, hj)) 



Z. 



(6) 



Distribution of partition function zeros of the ±J model on the Bethe lattice 



5 



limiting distribution satisfies the following self-consistent equation 



l[PH,p{hj)dhj. 



(9) 



PH,p{h)= / 6{h-H-J2''.^{h,)) 

I :i=i Ijr- 

The existence of a unique solution of this equation is shown rigorously in [13] for real 
values of the cavity field. Although this proof is not directly applicable to our case of 
complex cavity field, we observed numerically that equation (9) has a solution for p < 1. 
The special case of pure ferromagnet (p = 1) in a purely imaginary external field is 
discussed in detail in Appendix A. 



2.2. Definition of the Bethe lattice 

The Cayley tree is pecuhar because the number of surface sites is comparable to the 
total number of sites and hence the contribution from the surface cannot be neglected. 
This property is inconvenient for studying the bulk of the system. Thus the Bethe lattice 
is often considered instead. There are two different definitions of the Bethe lattice, and 
we have to clearly distinguish them [14]. 

(i) The interior part of the Cayley tree: A lattice consisting of the interior part of 
an infinite-size Cayley tree with uncorrelated boundary conditions. Alternatively, 
we can define this lattice as a finite Cayley tree, the boundary conditions of which 
are given by the convergent cavity field distribution Ph,is = lim;_^oo of the 
infinite-size Cayley tree. 

(ii) The regular random graph (RRG): A randomly generated graph under the 
constraint of a fixed connectivity c. Since there exist many cycles, we cannot exactly 
treat the finite system size. In the limit N ^ oo under appropriate conditions (and 
outside the spin-glass phase), however, the contribution coming from the loops is 
expected to become negligible and the problem can be solved exactly by the cavity 
method. 

Both definitions have advantages and disadvantages, and typically the phenomena 
that occur in one set-up have a correspondence in the other. The model (i) and its 
phase diagram have been discussed in detail in [13], while for (ii) and its relations with 
the replica theory we refer the reader to reference [12]. 

In this paper, we refer to model (i) (with typical boundary conditions which are not 
correlated among each other or with the couplings in the interior) and call it the Bethe 
lattice. Although definition (ii) is useful in that it has no ambiguity about the boundary 
conditions, it is incompatible with the formulation of the partition function zeros that 
uses equation (3). On the other hand, according to definition (i), we can easily define a 
finite-size system of the Bethe lattice and look at its zeros in the formalism of the cavity 
method. 
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2.3. The partition function zeros of the Bethe lattice 

Let us first consider the zeros of the partition function with respect to the external field 
H, the Lee- Yang zeros [1]. It is known that these zeros are related to the magnetization 
of the system. In general, the partition function of an Ising system in an external field 
is expressed as a two- variable polynomial, 

N Nt 

Z {H, P) = e^^^e^"'^^ ^) 6-2^^-^6-2^^^, (10) 

M=OE=0 

where N is the number of sites and Ni, is the number of interactions. We can write the 
above equation with fixed temperature as 

N 

Z{H) = Ce^^^n (^"'^'' - ^"'^""O' (11) 

i 

where ^ is a constant (to be ignored hereafter). Equation (11) means that the partition 
function is a polynomial of degree in the fugacity e'"^^^ , and there are roots on 
the complex H plane. Taking the logarithm and dividing by N, we find 

-/3f {H)^/3H + jj (PH'g {H') log (6-^^^ - 6-^^^') . (12) 

where gn is defined as the density of zeros on the complex H plane. The complex 
magnetization m{H) — —df/dH is expressed as 

^ (^) = 1 + // d'H'g {H') (13) 

We can express the density of zeros in terms of the magnetization by an infinitesimal 
closed fine integral, 

as is easily seen from the representation (13). Indeed the line integral picks up all poles 
within the circle, and all poles have the residue 1//3. 

For the Bethe lattice, due to the uniformity of the system, the disorder averaged 
complex magnetization is given by 

JJ d^h^'^^PP^p tanh {0"'^), (15) 

where the distribution of the central complex field P^p is calculated from the convergent 
distribution Ph,i} as 



Pt,ih^'')=JJ 



6{h^^^-H-J2^,{h,)) 
Inserting the equation (15) into equation (14), we obtain 



l[PH,^{h,)d''hj. (16) 



= 2^ 15s ^ / / '""^ (^^^'^ 
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J J\H-H*\<r 



[ I cfH'p'i! {-Kim 



\H~H'\<r 



^P'^> (7rz/2/3). 



(17) 



In the second line, we used the residue theorem under the condition that the radius r 
is sufficiently small. This result connects the density of zeros at H with the density of 
cavity fields at fixed external field H. If we iterate the cavity field population numerically 
for a given pair of complex values {H,I3), the distribution function of the cavity field 
yields the density of zeros in the if-plane through equation (17). 

The relation (17) can also be interpreted as follows. The whole partition function of 
the Bethe lattice Z^'^^ is formally the same as that of the Cayley tree (6), the difference 
only being whether or not we use the limiting distribution PH,p{h). This implies that 
the equation of zeros Z^'^^ = becomes identical toj 



Equations (17) and (18) show that the complex support of the zeros of the partition 
function can be obtained from the knowledge of the distribution of the field acting on 
the central spin. This exact relation represents an important advantage of considering 
the Bethe lattice as the interior part of the infinite-size Cayley tree. It is one of the 
main results of the present paper. 

We now turn our attention to the density of zeros as a function of the complex 
temperature. The condition (18) for the vanishing of the partition function is still valid 
for complex values of temperatures. In order to find the correct density one should follow 
steps analogous to those which lead from (11) to (17), by treating df/df3 instead of 
df/dH. We have however not followed this route here. Instead we contented ourselves 
with detecting the support of the zeros, i.e., the location on the temperature plane 
where g (/3) is non-zero. This can be performed by using the fact that the value of the 
density on the field plane at a complex temperature /3 — /3', g{H — if')|^^^,, should be 
proportional to the density on the temperature plane in a field H — H', g {f3 — ^')|^^^,. 
This means that the equality g {H = i^')l/3=/3' — ^ W) 9 = /^')\h=h' holds, where 
C(/3') is a normalization factor (to be dropped hereafter) §. This equality is sufficient 
to understand the phase diagram in terms of the locations of zeros, while a quantitative 
knowledge of the density of zeros on the temperature plane cannot be obtained due to 
the nontrivial normalization factor C(/3'). Anyway we do not discuss quantitatively the 
value of density on the temperature plane in the present work. 

I Note that equation (6) may seem to diverge when the factor cosh/Swj becomes 0, but this is not 
the case. The reason is that the condition cosh/3ztj = is always accompanied by Zj = and 
{cosh pUj)~^ Zj yields a finite value. 

§ We can accept this equality by considering fact that the zeros are located in the four-dimensional 
complex temperature-field space, and the two-dimensional distributions g with real H and g (H) 
with real /3 are just two-dimensional cross sections of the same distribution function in four dimensions. 



2 cosh (/3/i(^)) = 0^/3/i(") 




(18) 
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3. Numerical procedures 

Hereafter, we discuss the density of zeros g{H) on the H = 2(5 H plane since H always 
appears in this combination. A hat on a variable, like will mean the same quantity 
multiplied by 2/3. The location of the zeros on the complex temperature plane will also 
be analyzed. 

For the spin-glass model, the zeros on the complex field plane are not restricted 
to the imaginary field axis but a finite fraction spreads into the complex plane. The 
remaining fraction still lies on the imaginary field axis. The total density of zeros on 
the field plane thus naturally splits into two parts, 

g{H)^S{^H)g^{0)+g2{H), (19) 

where 9 — (— tt < ^ < tt) is the imaginary part of H. The one-dimensional measure 
gi is restricted to the imaginary axis and the two-dimensional measure g2 is continuously 
distributed on the complex plane. 

We base our analysis mostly on equation (17). However, it is not sufficient for 
our purpose because it only yields the two-dimensional part g2. In fact, the cavity 
field distribution always spreads in the two-dimensional /i-planc due to the initial 
conditions that we choose for the Bethe spin glass, as explained later. It is not 
possible to obtain the one-dimensional part gi together with the two-dimensional g2 
from equation (17) ||. However we can calculate gi by a different method using the fact 
that the one-dimensional density on the imaginary field axis corresponds to the jump 
of the real part of the complex magnetization as H crosses the imaginary axis. In other 
words, 2ngi (9) = ^ m{H = i9 + 0+) [1]. The analogy with electrostatics developed 
in reference [1] allows us to apply this formula even when g2 7^ 0. Therefore we can 
calculate gi from equation (15). 

Wc now obtain the distribution of P^^ (/i*-^-* = 2/3/i*^'^^) for complex field H and 

temperature T. Let us explain the algorithms to obtain P^^ (h^'^^) and estimate the 
density of zeros. We write the recursion relation (5) as 

Uj = tanh~^ ^tanh (^ J^^) tanh + "'^ j j j ' 

where the interaction Jjj connects the current site j with the next site i, and k labels 
the previous sites. The central cavity field is 

c 

= ^ + 2/3^Mj-. (21) 

Our numerical solution for the probabihty densities i^u {u) and P^^\U'^^) is based on 
the method of population dynamics. We represent the distribution functions tth (u) and 
(/?,('=)) in terms of a large number of variables {uj} and {^^'^^}, whose distributions 
are supposed to follow the respective probability distributions. The elements of these 

II For non-frustrated systems, we can obtain gi from the relation (17) (see Appendix A). 



Distribution of partition function zeros of the ±J model on the Bethe lattice 



9 



sets are updated by randomly choosing Ui and Jij and following equations (20) and (21). 
The initial condition of the cavity bias, which correspond to the boundary condition of 
the Cayley tree, is taken to be uj = limH-^ooUj = Jij (g M). These initial fields are 
uncorrelated random variables, which introduce frustration into the system. Note that 
in general the value of Uj has both real and imaginary parts because the value of the 
next generation is calculated from equations (20) and (21) with complex j3 and H. 
The distribution of h^'^'' is calculated from the Monte Carlo (dynamical) average of the 
population in the two dimensional complex plane, after the dynamics has converged to 
a limiting distribution. 

The actual steps of the algorithm are as follows; 

(i) Set the initial probability distribution of Uj as Uj — Jij. 

(ii) Update the sets of {uj} and {/i^"^^} with the recursion relations (20) and (21) by 
randomly choosing Uk and Uj out of the set {uj} until P^^ {h^'^^) converges. 

(iii) Estimate the density of zeros by 

g^{e,T) = m{H = ie + 0+) 

= \im^ J I d^U'^'^P^f' (h^^^^ tanh (h^'^ 12 

g^{H,T)^pf{U^^ ^m) (23) 

for given H and T. The one-dimensional density gi represents a density of zeros 
under a pure-imaginary field while the two-dimensional density g2 is defined on the 
whole complex planes. 

For simphcity only the Bethe lattice with connectivity c — ?> has been studied. We 
have chosen A/pop = 10^ representative points in the population dynamics and have 
performed at least 5000A^pop cavity iterations until the population converges. The 
data were collected from the average of additional 5000A^pop iterations after the initial 
SOOOA^pop (or more) iterations. These conditions have been used throughout this work. 



(22) 

H=ie+Hjt 



4. Distribution of zeros for the Bethe spin glass 
4.1. Zeros on the complex field plane 

First, we show the density of zeros on the complex H plane for real temperature. Figure 
2 is the distributions of zeros on the complex H plane with probability p = 0.5 at 
temperature T = 1.43 > Tsg = l/tanh"^ (l/Vc^) (left) and T = 0.5 < Tsg (right). 
The complex plane has been spht into cells by dividing the real axis from 3? (H) = to 
12 with an increment 0.25 and the imaginary axis from $5 (H) = 0.02 to 7r/2 with an 
increment of 0.02. The density outside this range is omitted, since this is sufficient to 
see zeros near the real axis which are essential for critical phenomena. Both gi and (72 
are plotted in the same figure and coloured in a logarithmic scale; a black dot shows a 
very high density. 




5 10 5 10 

Figure 2. Distribution of zeros on the complex H plane with p — 0.5 at T — 1.43 
(paramagnet, left panel) and T = 0.5 (spin-glass, right panel). Densities are coloured 
in a logarithmic scale. On the right panel, zeros with finite real part touch the real 
axis, whereas gi on the imaginary axis vanishes to the numerical precision. 



The left panel of figure 2 (T = 1.43) is for the paramagnetic phase. Neither gi nor 
g2 has a finite value on the real axis and there is thus no phase transition as a function 
of the real field. The right panel of figure 2 is in the spin-glass phase (T = 0.5), where 
zeros reach the real axis at and away from the origin. Thus, there is a phase transition 
at some H = Hc{> 0) G M on the Bethe lattice, where is the point where the density 
g2 changes its value from zero to non-zero along the real axis. Zeros touch the real 
axis all the interval between H = and H = He- This suggests that the free energy 
is non-analytic as a function of H below He- The one-dimensional density gi vanishes 
to the numerical precision for both temperatures in figure 2. The two-dimensional part 
g2 has finite values on the imaginary axis. It is likely that the edge of g2 (the point 
where (72 7^ on the imaginary axis and is closest to the origin) touches the origin as the 
temperature is decreased to T = Tsg- At lower T the density of the zeros approaches 
the real axis also away from the imaginary axis, up to a point Hc{T) on the real axis. 
This defines a line of spin-glass transitions in the H-T plane, as discussed below. 

In the ferromagnetic phase (the right part of figure 3), on the other hand, only gi 
touches the origin but g2 does not. In marked contrast to the left panel of figure 2, the 
one-dimensional density has a finite value on the imaginary axis away from the origin on 
the left panel of figure 3. This is an important difference of ferromagnetic and spin-glass 
transitions seen from the distribution of Lee- Yang zeros. More details on the transition 
points will be discussed in later sections. 
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5 10 5 10 

Figure 3. Distribution of zeros on the complex H plane with p = 0.9 at T = 1.5 (left: 
para) and T = 0.5 (right: ferro). Only the one-dimensional zeros on the imaginary 
axis, gi (thin lines), touch the origin at low temperature. 



4-2. Zeros on the complex temperature plane 

Using complex values of the temperature in equation (20), we can also obtain information 
on the support of Fisher zeros, using equation (23). This information is sufficient for 
our purpose of detecting phase transitions from the point of view of zeros. 

Figure 4 shows the distribution of zeros on the temperature plane at p = 0.5 with 
H = (left) and H = 0.5 (right). The (two-dimensional) zeros also touch the real 
temperature axis below the spin-glass transition temperature. This confirms that the 
spin-glass transition differs from the ordinary phase transition where zeros touch the 
real axis only at the transition temperature. This may be interpreted as the system 
staying critical at all temperatures below the transition temperature. It may also be 
taken as a signature of temperature chaos, i.e. the instability of the randomly frozen 
spin configurations in the spin-glass phase with respect to arbitrarily small temperature 
changes. In accordance with the right panel of figure 2, the right panel of figure 4 
shows that the spin-glass phase persists under a weak field {H = 0.5). The critical 
temperatures appearing in figure 4 are consistent with the phase diagram shown in the 
next subsection. 

The distributions at p close to 1 are also interesting. Figure 5 shows the density at 
p = 0.9 with H = (left) and H = 10~'^i (right). The transition temperature between 
paramagnetic and ferromagnetic phases is known to be Tc ~ 1.36. It seems that the zeros 
distribute on an almost one- dimensional curve near which most likely touches the 
real axis at T = Tc. Since our method of equation (23) calculates the two-dimensional 
density of zeros, a one-dimensional density is hard to identify with a high precision. On 
the right panel of figure 5, a weak field H = lO^H is added, to test for a spontaneous 
magnetization: zeros have finite densities along the real axis below T^, since gi is finite 
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Figure 4. Distribution of zeros on the complex T plane with p = 0.5 and H = 
(left) and H = 0.5 (right). Zeros touch the real axis below the critical temperature 
TsG — 1-13 (left) and 0.5 (right). The apparent absence of zeros near the origin may 
be due to numerical rounding errors of tanh/3. 



under a pure imaginary field in the ferromagnetic pfiase. 

On tlie otlier liand, tlie two-dimensional distribution g2 touches the real axis again 
in the low temperature region. Below Tsg — 0.29, the zeros approach the real axis on 
the whole interval < T < Tsg, similarly as in the spin-glass phase shown in figure 
4. Indeed, this second critical temperature corresponds to the spin-glass transition, as 
shown below. 

4-3. Phase diagram 

Based on the above observations, we investigate the phase diagram on the p-T and 
T-H planes with real T and H. Graphical representations shown so far suggest that 
there are two types of transitions where one- and two-dimensional distributions of zeros 
reach the real axis separately. In order to determine the two transition points, we 
added a very small imaginary field H = 10~H. Since zeros on the imaginary axis in 
the complex H plane reach the real axis before those away from the imaginary axis as 
the temperature is decreased, we can restrict to the imaginary H axis and ask for the 
highest temperature for which a non-zero density, gi and g2 exists at the origin of the 
imaginary axis {H = 10~^i). In contrast, it is difficult to determine the transition points 
caused by one-dimensional zeros on the complex T plane in the absence of a field as 
shown in the left panel of figure 5. The transition temperature ~ 1.36 derived from 
the first method of the complex H = 10~^i is close to the point in the left panel of 
figure 5 where the almost one-dimensional distribution is likely to touch the real axis. 
Therefore we mainly used the method of complex field to identify the phase boundary 
(the left panel of figure 6). 




0.29 7^=1.36 0.29 7^=1.36 

Figure 5. Distribution of zeros on the complex T plane for p = 0.9 and H = (left) 
and H = lO^^'i (right). The zeros approach the real axis at T = Tc, where Tc ~ 1.36 
is the ferromagnetic critical temperature. In addition, zeros reach the real axis also in 
the low temperature region. The second critical temperature corresponds to the spin- 
glass transition Tsg shown in figure 6. On the right panel, a weak pure-imaginary 
field immediately makes apparent the presence of ferromagnetic order below a critical 
temperature T^. a one-dimensional distribution of zeros lies on the real axis in the 
right panel. 

Our result agrees well with the exact phase boundary [15, 16, 17, 18], Tc = 
l/tanh~^ [1/ (4p - 2)] (between para and ferro phases) for p > Pc = (2 + a/2)/4 and 
Tsg = 1/tanh"^ [1/^2] (between para and spin-glass phases) for p < p^. It is expected 
that the one-dimensional distribution gi determines the ferromagnetic temperature Tc 
and the two-dimensional distribution g2 determines the spin-glass transition temperature 
Tsg- We discuss this hypothesis in detail below. 

4.3.1. Behaviour of the two-dimensional distribution of zeros 

Below the line drawn in circles in figure 6, the ordered phase is expected to be stable 
under an external field as shown in figure 2. This line is also determined from the zeros 
on the temperature plane. We set the temperature as T = Tr + lO^^i with H = 0, 
where Tr is real, and we assume the highest Tr having a non-zero density as the critical 
temperature (triangles in figure 6). The two results agree well with each other. 

The boundary between the ferromagnetic and spin-glass phases is harder to 
determine. As one sees in figure 6 (the left panel), the temperature where the one- 
dimensional distribution gi ceases to touch the real axis is lower than the temperature 
where the two-dimensional distribution (72 starts touching the real axis (for biases 
p ~ 0.85 to 0.92). This indicates that there is a phase where the system is still 
spontaneously magnetized, but also glassy (see also [13]). It is shown in Appendix 
B that the spin-glass susceptibility xsg diverges along a line which coincides with the 
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Figure 6. Phase diagrams on the p-T plane from the density of zeros. The solid lines 
are analytical results as given as Tc = 1/ tanh^^ [1/ (4p — 2)] (para-ferro boundary) 
and = l/tanh~"^ [l/'v/S] (para-spin glass boundary) [15, 16, 17, 18]. Left: Circles 
and triangles are determined by the two-dimensional density of zeros 172 on the complex 
field plane and temperature plane, respectively, which corresponds to spin-glass critical 
temperatures. Squares are phase boundaries calculated by the one-dimensional density 
of zeros 171, which indicates the onset of the ferromagnetic order. Immediately below 
triangles there is a phase which is both magnetized and glassy. The dotted line denotes 
the Nishimori line [19], on which the multicritical point (where three phases merge) 
is located. Right: Blow up of the field-induced critical temperatures on the left panel 
(marked in circles and triangles) between p = 0.84 and p — 0.93. 



touching points of g2, marked in circles and triangles in figure 6 . This implies that 
our method using the approach of the two-dimensional distribution (72 to the real axis 
correctly detects the onset of the spin-glass phase. Indeed the proximity of zeros to the 
real axis at if = implies the divergence of some (higher order) susceptibility. 

We next show the phase diagram on the T-H plane with fixed p (figure 7). We 
added a complex external filed as H = + 10~^i and we assumed that the critical 
value of He is the maximum value of HR/2f3 = H having a non-zero (two-dimensional) 
density g2- The top panels of figure 7 are for the temperature dependence of the critical 
line for various p. For all biases p, smoothly rises from if = at T = Tc {H = 0) for 
p < Pc= (2 -f a/2) /4 ^ 0.854, which is for the multicritical point on the Nishimori line 
[19] where the three phases merge. At low temperature, it seems that He approaches 
a finite value at T = 0. The bottom panel is for p = 0.5, and the points well fit 
to He oc (Tc — T)^'^. The exponent 1.5 is the same as in the AT line [20] of the 
Sherrington-Kirkpatrick model [21]. The existence of the AT line on the Bethe lattice 
(regular random graph) was predicted in reference [15] where the AT line should behave 
like (Tc — T)^^"^ in the vicinity of T = Tc, and was checked by the population method in 
references [22, 23]. 

The agreement of the exponent 1.5 with the behaviour of the AT line might be 
interpreted as an indication that the phase below the phase boundary is the spin-glass 
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Figure 7. Spin-glass transition 
lines in the T-H plane (corre- 
sponding to the AT line in the 
SK model). The bottom panel 
is for p = 0.50. The boundary 
at p = 0.5 rises from T — Tc 
as He cx (Tc-T)^-^. The expo- 
nent 1.5 agrees with that of the 
AT line for the SK model. Pc = 
(2 -I- V2) /4 ~ 0.854 is for the mul- 
ticritical point. 



phase with rephca symmetry breaking (RSB) [24, 25, 26]. Indeed, it is suggested that 
the spin-glass phase in the regular random graph has the full RSB [27, 28]. However, we 
should be careful because we have not directly seen the instability of a replica-symmetric 
solution in the Bethe lattice of our definition. The analysis of reference [13] shows that 
in the spin-glass phase a system of two replicas on the Cayley tree exhibits a diverging 
susceptibility with respect to an infinitesimal repulsion between the replicas. However, 
at the same time the cavity field distribution remains essentially identical to the replica 
symmetric approximation for the regular random graph; model (ii). Further studies of 
this phase from different viewpoints may be necessary to fully clarify its nature. 

4.4- Griffiths singularity 

The density of zeros on the imaginary field axis is closely related to the Griffiths 
singularity in the diluted ferromagnet [3]. The same is expected in spin glasses [8]. 
The Griffiths singularity is expected to manifest itself, if it is present, in the form of 
an essential singularity of the density of zeros upon approaching the origin along the 
imaginary field axis [4, 5, 6, 8]. If such a tail is present (which is very difficult to detect 
numerically), its touching of the real axis would indicate the onset of a Griffith phase. 
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Figure 8. One- and two-dimensional density of zeros on the imaginary field axis for 
high p at T = 1.5. The top panels are for the ferromagnetic phase and the bottom for 
the paramagnetic phase. The contribution of the one-dimensional density gi decreases 
with decreasing p whereas 52 increases. 



but not the phase transition to a ferromagnetically ordered or glassy state. In the 
presence of a Griffith phase, our criterion for the detection of phase transitions should 
thus strictly speaking be refined to the condition that a "substantial density of zeros" 
(which grows at least like a power law with away from the real axis) touches the 
real axis. 

The above discussion suggests to study in more detail the form of the one- 
dimensional density of zeros on the imaginary axis as an indicator of possible Griffiths 
singularity for the Bethe spin glass. 

For the Bethe spin glass our numerics does not show detectable signs of a Griffiths 
singularity above the spin-glass phase, our estimated gi being zero within numerical 
precision in the paramagnetic phase above the spin-glass phase, as shown in figure 2. 
Above the ferromagnetic phase, it seems that gi is finite above some Oq (> 0), while 
there seems to be a jump between gi = {9 < 9o) and > > ^o)- 

Figure 8 shows the one- and two-dimensional densities of zeros on the imaginary 
axis at fixed T = 1.5 with various p. The top panels correspond to biases p which are 
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in the ferromagnetic phase while the bottom panels correspond to lower biases, in the 
paramagnetic phase. At high p, the gi part is dominant because the line integrals 



for p = 0.99 and 0.925 are nearly equal to 1, where Wi + W2 = 1 with W2 = 



decreases with decreasing p, and gi vanishes for large 6. In the left-bottom panel of 
figure 8, gi is seen to decrease when g2 becomes finite around 9 = 1.0. This behaviour 
suggests that one-dimensional zeros are buried under large cloud of two-dimensional 
zeros far away from the real axis. Finally, for p < Pc the one-dimensional gi seems to 
vanish altogether above the spin- glass phase as shown in figure 2. 

It is difficult to draw a definite conclusion about Griffiths singularities due to the 
limit to detect very small but non-vanishing values of gi near the origin. However, our 
numerics points toward the possible absence of a Griffiths phase in the spin glass on 
the Bethe lattice, contrary to what happens for finite-dimensional systems and diluted 
ferromagnets. A more detailed study will be required to settle this question. 

5. Conclusion 

We have investigated the distribution of the partition function zeros for the ±J spin- 
glass model on the Bethe lattice. An important relation to connect the density of 
zeros and the density of cavity fields for complex field and temperature is found, which 
enables us to treat the zeros in the limit of infinite system size. The densities are split 
into one- and two-dimensional parts, where the one-dimensional density is defined on 
the imaginary axis of the complex field plane and the two-dimensional density spreads 
over the complex plane. We investigated the phase diagram of the Bethe spin glass 
by estimating the transition points where the one- or two-dimensional densities begin 
to have finite values in the immediate vicinity of the real axis. Our results agree well 
with the analytically exact phase diagram. The one-dimensional density determines the 
boundary between the paramagnetic and ferromagnetic phases since the one- dimensional 
density is directly related to the spontaneous magnetization. The two-dimensional 
density determines the boundary of the spin-glass phase on the p-T and H-T planes. 
This phase boundary corresponds to the line where the spin-glass susceptibility diverges 
as evaluated in Appendix B. Below the spin-glass transition point, the two-dimensional 
density of zeros continuously touches the real field and temperature axes, which may 
be related to chaotic behaviour of the system as a function of the field and temperature 
in the spin-glass phase. We have shown, by looking at the locations of the partition 
function zeros, that the system has a non-analytic free energy for all T below the spin- 
glass temperature Tc. 

We have not observed any evidence for the existence of a Griffiths phase in our 
numerics. If such a phase is indeed absent, it implies that the Bethe lattice behaves 




(24) 



JJ dHg2 {H )dH should hold due to the normalization. On the other hand, Wi rapidly 
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distinctly from finite-dimensional spin glasses. Further careful studies are required to 
understand the origin of differences between Bethe and Bravais lattices. 
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Appendix A. Zeros of the pure ferromagnetic system 

In order to check our method itself and the precision of the numerical analyses, we 
estimated the density of zeros for the pure ferromagnetic system. The zeros of a 
ferromagnetic system is located only on the imaginary axis, and thus g2{H) = 0. In this 
case the density of zeros on the imaginary axis, gi{H), can be found by both methods 
which we used to calculate gi and g2 in the case p < 1 (equations (17) and (22)). We 
here compare the exact density of zeros and two algorithms using equations (17) and 
(22). Figure Al is for comparison of numerical calculations and the exact solution. Our 
algorithms agree perfectly with the exact solution. The details of the analyses are as 
follows. 

(i) The exact density of zeros (the solid line in figure Al): The density of zeros on 
the imaginary field axis gi {9) is obtained as the analytic continuation of the real 
m{2/3H) from real positive values oi H — 2/3H to purely imaginary H — iO, 

2'ngi {e) = ^m {H := iO). (A.l) 

Using the fixed point hf {H,(3) of equation (5) with J^j = 1 and c = 3, we rewrite 
equation (A.l) as 



2TTgi (9) = 3? tanh 



(A.2) 



H:=ie 



From equation (5) the fixed point hf {H,^) for the pure ferromagnetic system in 
real field is found as 

hf{H, (3) ^2logx', (A.3) 

where x' is the solutions of the following equation, 

e^x' - e^+^^x^ + e^^x - 1 - 0. (A.4) 

(ii) The relation between density of zeros and cavity fields (squares in figure Al): The 
algorithm shown as equation (17) connects the density of cavity fields and density 
of zeros. Since the zeros for the ferromagnetic system lie only on the imaginary 
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field axis, the density gi (9) is estimated from the distribution of cavity field in a 
purely-imaginary field. We thus consider the free boundary condition in order to 
make the cavity field pure imaginary. 

For the cavity field in a pure imaginary field 6, the recursion relations (20) and (21) 
are rewritten as 

Qui = 2tan-i [tanh/3tan(^/2 + (c- l)$>%/2)] ^Ui = 0, (A.5) 

c> /iW = e + c^Ui 3? h^''^ = 0, (A.6) 

where it = 2f3u and h^'^'^ = 2f3h^'^^ and the initial u is set as -u = 0. We directly 
estimate the one-dimensional density of the population on the imaginary field axis 
by using the relation, 

g,{e) = Pi^\'^k = n). (A.7) 

It can be shown in full generality that this procedure yields a density gi{H) 
which correctly describes the jump of the real part of the magnetization across 
the imaginary H axis. For the spin-glass model {p < 1), this algorithm is used only 
in the estimation of g2 part (equation (23)). Since the boundary condition of the 
Bethe spin-glass is fixed in order to introduce frustration to the system and the 
cavity fields are complex, this method is not applicable to the estimation of gi for 
the spin-glass model. 

(iii) The real part of complex magnetization (circles in figure Al): Complex 
magnetization in the vicinity of the imaginary field axis is estimated from 



m (9) — lim tanh 



(A.8) 

H=ie+HR 



where h'j^ is the (numerically) fixed point of the complex cavity field. The one- 
dimensional density on the imaginary axis is also calculated from the real part of 
the complex magnetization as gi (9) = {9) /2'k. This method is used to calculate 
gi for the Bethe spin-glass (equation (22)). 

Note that the imaginary cavity field (equation (A.5)) docs not converge for the pure 
ferromagnetic system when the density of zeros has a positive value for a given H (the 
left panel of figure A2). However we take the average over a large number of the iterating 
steps to calculate the density of cavity field, which naturally performs the sampling of 
bulk properties of the Cayley tree. If the return map of the recursion relation (5) has a 
fixed point, on the other hand, the cavity field population is delta-peaked at that fixed 
point. One finds that the density of zeros vanishes under this condition (see the right 
panel of figure A2). 

The critical value of the imaginary field is given by the condition that the two 
curves in figure A2 touch: 

2 (c - 1) tan"^ tanh p tan (9o/2 + h/2^ (A.9) 

and 

i [2(c-l)tan-^ [tanh^tan(^o/2 + /i/2)]] =1. (A.IO) 
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Figure Al. The density of 
zeros on the imaginary axis for 
the pure ferromagnetic system 
with c — 3. The numerical 
estimations are consistent with the 
exact calculation of equation (A. 2). 
The edge of zeros corresponds to 6*0 
in the text. 
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Figure A2. The return map of equation (21) for the pure ferromagnetic system in 



4 tan 



tanh j3 tan 



((e + S/i) /2 



does 



pure imaginary field. Left: If the curve y 

not intersect y — the imaginary cavity field does not converge; the density of 
zeros has a finite value under this condition. Right: If the two equations intersect, the 
imaginary cavity field recursion has fixed point; thus the density has no value. 
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We find that the critical imaginary field is determined as 



-1 



^0 = 2 cos 
-2 (c - 1) tan-i 



a/(c — 1 — tanh /3) sinh /3 cosh /3 



(1 - (c- l)tanh/3) tanh/3 
c — 1 — tanh (3 



(A.ll) 



The critical temperature is given from this equation at 6'o = as /3 = tanh^ [1/ (c — 1)]. 
This result is of course consistent with the well-known critical temperature. Figure A3 
shows the temperature dependence of 9o for c = 2, 3, ■ ■ ■ 9. 




Figure A3. The temperature 
dependence of 9q for connectivities 
c = 2, 3, • • • 9 from left to right. At 
^0 = 0, the temperature is given by 
T = l/tanh-^[l/(c-l)]. 



Appendix B. The critical condition based on the spin-glass susceptibility 



It is quite possible that the divergence of the spin-glass susceptibility XsG, which in 
the SK model is identified with the AT instability [14, 29], characterizes the instability 
signaled by the zeros of the partition function approaching the real axis. Here we show 
that this is indeed the case for the present Bethe lattice. 
The spin-glass susceptibility is defined as 



XSG 



N ^ 



d{S,) 
dh, 



d{So) 
dh^ 



fB.ll 



To derive the second identity, we assumed uniformity of the Bethe lattice and selected 
the central spin as i. In a cycle-free graph, an arbitrary pair of sites are connected 
by a single path. Let us assign site indexes from the origin of the graph to a site of 
distance G along the path as g = 1,2, ... ,G. For a fixed set of couplings and boundary 
fields, the chain rule of the derivative shows that 



d {So) d {So) dho duo dha 



G 



dh, 



G 



du 



G 



d {So) dho jj dug 



-1 dhn 



dhr, dur 



9=1 



dh„ du. 



dh^ duQ dh^ 

_ d {So) -pr dUg^i 

duo f}^ dUg ' 

as hg is the cavity field on site g and depends linearly on Ug as hg 
rg represents the sum of the cavity biases from the other branches that fiow into site g. 



(B.2) 



H + Ug + rg, where 



Distribution of partition function zeros of the ±J model on the Bethe lattice 



22 



In the hmit G — >■ oo, the relevant factor in equation (B.2) is only 11^=1 i^Ug-i/dug) . 
The spin-glass susceptibihty is hence evaluated as 

2" 



XSG = Yl 



G-1 



G=l 



d{So) 



dh 



G 



oo 



J G=l 



■ G 

n 



du„ 



,(B.3) 



where the factor c(c — l)*^ denotes the number of sites of distance G from the central 
site 0. The divergence condition of xsg is given by 



log(c - 1) + hm — log 

G^oo G 



■ G 

n 



0. 



(B.4) 



This yields the condition of the spin- glass transition of the Bethe lattice and RRG. 
In order to estimate the divergence points of the spin-glass susceptibility, we 



numerically implement the calculation of the factor Yl„=i(dug-i/dug) . This factor 

L ^ J J 

can also be written as [{duo/duG)^]j and this latter form is more tractable at finite 
temperatures. 

The numerical evaluation of the factor duo/duc is straightforward, 
duo uo {ug + Amg) - uq {ug) 



"G ' ^^'^^ 

The procedure to evaluate this equation is as follows [28, 30, 31]. Wc arrange two 
replicas of an identical population {ui}^^l^ expressing the convergent (real) cavity-bias 
distribution tih,^ {u), which is related to the convergent cavity- field distribution PH,j3{h) 
given in equation (9) as 



t^hA'^) = J dhPnA^) 



8 [u — — tanh ^ {tanh ^ J tanh ^h} 



(B.6) 



In addition, we introduce a uniform perturbation (Am = 10""^) into only one of the two 
rephcas and then observe the square average of the variation, (1/A^pop) Yl!i=\'{^ik^G + 
AiiG) ~ 'Ui{uG)Y 1 s-fter a certain number of the cavity updates. 

In particular, we update two populations by SOOOA^pop iterations with the same set 
of Jy. A critical line of the divergence of the spin-glass susceptibility is determined by 
whether the square average is numerically zero or much larger than the perturbation. 
The result is shown in figure Bl where the spin-glass susceptibility diverges below 
this line. At zero temperature, the critical probability is calculated as p,. = 11/12 in 
references [27, 32], which is reproduced by our numerical calculation at zero temperature 
aspc = 0.91665(5). 

This result agrees well with the phase boundary drawn in figure 6. Thus it is 
suggested that our phase boundary estimated by the two-dimensional zeros corresponds 
to the spin-glass transitions. 
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